Cooling and aggregation in wet granulates 
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Wet granular materials are characterized by a defined bond energy in their particle interaction 
such that breaking a bond implies an irreversible loss of a fixed amount of energy. Associated with 
the bond energy is a nonequilibrium transition, setting in as the granular temperature falls below 
the bond energy. The subsequent aggregation of particles into clusters is shown to be a self-similar 
growth process with a cluster size distribution that obeys scaling. In the early phase of aggregation 
the clusters are fractals with Di = 2, for later times we observe gelation. We use simple scaling 
arguments to derive the temperature decay in the early and late stages of cooling and verify our 
results with event-driven simulations. 
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Granular systems have encountered strongly increasing 
interest within recent years [l|, Q as they are both close 
to industrial applications and relevant as model systems 
for collective systems far from thermal equilibrium. Most 
research has focused on dry granulates resulting in rich 
non-equilibrium phenomena The dramatic change of 
mechanical properties due to some wetting liquid addi- 
tive is apparent, when we compare the fluid-like state of 
dry sand flowing through an hourglass with the shapeable 
plastic state of wet sand, which is suitable for molding 
a sandcastle. This change in bulk properties can be at- 
tributed to the differences in the underlying particle in- 
teractions [1, [si]. Whereas in collisions of dry particles, a 
certain fraction of the initial kinetic energy is dissipated 
into the atomic degrees of freedom of the particles, in the 
wet case, the interaction is mainly due to the interfacial 
forces exerted by liquid capillary bridges which form be- 
tween adjacent particles and dissipate energy upon rup- 
ture. 

More recently, the dynamics of wetCTanular media has 
been addressed in several studies d, [a0, IMJl, , focus- 
ing on nonequilibrium phase transitions 10'|, agglomer- 
ation 0| , shear flow Q and cooling in one dimension 
0. A particularly important aspect of free cooling in 
cohesive gases is the structure of the emerging clusters. 
It pertains to the formation of dust fllaments and the 
microscopic mechanisms of cloud formation as well as 
to the size distribution and impact probability of plan- 
etesimals in accretion discs. Structure formation in wet 
granulates during free cooling has hardly been studied 
yet and is the central theme of our paper. We suggest a 
very simple model for the interaction of two wet grains, 
which only takes into account the essential features of a 
capillary bridge: hysteresis and dissipation with a well 
defined energy loss. Cooling is controlled by the proba- 
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bility for a bridge to rupture and hence logarithmically 
slow in the long time limit, when a percolating structure 
has been formed. For smaller times the structure is char- 
acterized by coexisting fractal clusters of all sizes, whose 
size distribution is shown to scale. 

Model — Consider a gas of N hard spheres of diameter 
d in a volume V — L^. Each particle is covered with 
a thin liquid film. When two particles collide, a liquid 
bridge forms between the two spheres and induces an 
attractive force by virtue of the surface tension of the 
liquid [1, @ . As the spheres withdraw from each other 
and their separation s increases, the bridge continues to 
exert a capillary force, F(s), up to some critical distance, 
dc > d, where the bridge ruptures and liberates the par- 
ticles. Thus there is a hysteretic interaction between the 
grains. Each rupture of a liquid bridge gives rise to dis- 
sipation of a fixed amount of energy, /S.E = J^" F{s)ds. 
Note that during each collision the total momentum of 
the impact partners is conserved. If the relative velocity 
of the impacting particles is insufficient to rupture the 
liquid bridge, v < Va-it ~ \/'2AE/ii (where fj, is the re- 
duced mass), they form a bounded state and continue 
their motion with the center of mass velocity. 

We have simulated the above model with an event 
driven algorithm, assuming the force between two col- 
liding particles to be zero, except when they touch - 
such that their relative distance is equal to d - or re- 
cede beyond the critical distance dc for bridge rupture. 
When the particles touch, they are elastically refiected 
(complete restitution). When the colliding particle move 
apart and cross the critical distance dc, the energy AE 
is taken out of the motion of the impact partners - pro- 
vided V > t;ciit- Otherwise the particles are elastically 
reflected and remain in the bounded state. In a realistic 
model of wetted grains (thin film), the capillary film is 
infinitcsimally thin, so that it is justified to assume that 
a bridge is only formed, if the particles actually touch. 
This implies that energy is removed at dc only, if the 
grains have touched each other before. In an alternative 
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model (thick film) one assumes that a bridge is formed, if 
the particles approach closer than dc- The difference be- 
tween the two models is insignificant for our results and 
is discussed in more detail in [11] . We use dimensionless 
units such that AE — 1, particle mass m — 1 and particle 
diameter d = 4. The bond-breaking distance is chosen as 
dc = 1.07c? and volume fraction is varied from (f> sa 0.06% 
up to 15.6%. We use periodic boundary conditions in the 
X- and y-direction and hard walls in ^-direction. 

Cooling — We define the granular temperature T = 
3^ "^'^f investigate its decay in time from a 

given initial value Tq 2> AE. In each collision an amount 
of energy AE is lost with probability Pbb, the probabil- 
ity to break a capillary bridge. Hence the average loss of 
energy per time is given by 

l'§ = -l-f.on-AE.nb. (1) 

where /coii denotes the collision frequency. In the early 
phase of cooling T ^ AE, so that every collision also 
leads to a rupture of the created capillary bond and Pbt ~ 
1. The collision frequency, /coii = 'ig[d)any/T/'Km is 
well established for a dilute system, with the particle 
density n = y, the scattering cross section a and the 
pair correlation function at contact g{d). For the latter 
we use the Carnahan-Starling approximation (l^ . The 
simplified cooling equation ^ ^ —Vt is easily solved by 

with a characteristic time scale 

° " 2gid)anAE ' ^ > 

In Fig. [1] the evolution of the granular temperature T 
from the simulation (dots) is compared to the analytic re- 
sult (full lines) for different volume fractions (j). The main 
dependence on </> is contained in the timescale io- Hence 
the different curves can be approximately collapsed onto 
a single curve, if T is plotted versus scaled time t/to (see 
inset). 

The agreement between theory and simulation is quite 
good up to to, when the assumption contained in eq. ([T]) 
that every collision causes an energy loss AE breaks 
down. Yet the timescale tg has a clear physical rele- 
vance. At <o the temperature is comparable to the bond- 
breaking energy AE and for t > to persistent clusters 
will form with important consequences for the cooling 
dynamics. Thus the timescale to separates two regimes 
with qualitatively different structural and dynamic prop- 
erties. 

Aggregation — For t > to the system starts to form ag- 
gregates, which grow in a self-similar process, reminis- 
cent of ballistic cluster-cluster aggregation [l^ . For even 
longer times we observe a spanning or percolating clus- 
ter for all finite densities and ultimately all particles and 
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FIG. 1; (color online) Decay of the granular tem- 
perature T for volume fractions (from left to right) 
(j> = 15.6%, 7.81%, 3.90%, 1.95%, 0.98%, 0.49%, 0.24%, 0.12%, 
0.061%;number of particles A'^ = 262144 is fixed; comparison 
of theory (full line) and simulation (•) inset: T versus scaled 
time t/to- 




FIG. 2: Snapshot of the system with <j> = 0.98%, N = 262144 
and T ~ 0.06Ai5 taken at time t « 36to; only particles in the 
spanning cluster are shown. A rotating version of this figure 
gives a better impression of the three-dimensional structure 
(download at [ll]). 

clusters have merged into a single cluster. In Fig. [5] we 
show such a cluster which has attracted 99.7% of all par- 
ticles. 

The structure of the aggregates in the early stage of ag- 
gregation, sometimes called flocculation regime, is char- 
acterized by the fractal dimension D{. It is easily deter- 
mined from the radius of gyration, which is expected to 
scale with the cluster mass m like rg ^ m^/-^', defining 
Df. In Fig. [3] we show the radius of gyration for a system 
of 262144 particles at volume fraction (j> — 1.95%. Sev- 
eral snapshots of the ensemble of growing clusters have 
been taken at times to < t < 4to (nonpercolating regime) . 
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The data scale well according to the above relation, some 
scatter is observed for the largest masses, corresponding 
to times close to the percolation transition. Since our 
clusters are not completely static, we plot the fractal 
dimension for each snapshot separately in the inset of 
Fig. [3l The variation is very small so that we conclude 
that Df does not depend on time. 



50 

20 
10 



10 



2 






D, 1 





i *.\ ^ 


1 2 3 4 
f/'n J 











10- 10-^ 
mass m 



10" 



105 



10" 



0.01 



0.1 



10 



mim 



FIG. 4: (color online) Rescaled cluster size distribution 
f{m/m) — Nm{t) ■ versus scaled cluster mass m/m for 
(/> = 1.96%; color coding refers to different snapshots. 



I log 01, can we expect to find a nontrivial Df such that 
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FIG. 3: (color online) Radius of gyration as a function of 
cluster size for a system of 262144 particles at volume fraction 
(j> = 1.95%; different colors correspond to simulation times 
between to (yellow) and 4to (black) ; the solid line has a slope 
of 1/2; inset: fractal dimension for each snapshot separately. 

All information about the connectivity of the clusters 
is contained in the cluster size distribution, Nm{t), the 
number of clusters of size m at time t. For aggregating 
systems the mass distribution is expected to evolve to- 
wards a self-preserving scaling form [e.g.flH. independent 
of the initial distribution: 

7V,„(t)=m-V(m/m(i)), (4) 

where the time dependence is only contained in the mean 
cluster mass fh(t) = y^° ^_-^ m?Nm{t)/N . Mass conserva- 
tion requires 9 = 2 [l5| | . We plot in Fig. 0] the scaling 
function f{m/m) = Nra{t)m? for the data from eight 
snapshots taken at times t between and 2<o- The data 
are seen to scale well in this regime, in which the mean 
cluster mass increases roughly by a factor of 30. 

The asymptotic cluster, which contains all the parti- 
cles, is compact on large scales. However, its structure 
at smaller length scale may well be fractal. To inves- 
tigate the Hausdorff dimension of the largest cluster at 
intermediate length scales, we use the box counting algo- 
rithm [l7j and determine the number of boxes iVbox 
of edge length Lbox, needed to cover the whole cluster 
[r§ |. On length scales much smaller than the particle 
diameter, Lbox d, the system obviously behaves three- 
dimensionally and A^box = 4>{L / L\,o^)'^ ■ Since our system 
is finite and contains a system-spanning cluster, the scal- 
ing behavior on large length scales should also be three- 
dimensional: A^box = (L/Lbox)'^ for Lbox ~ L. Only in 
the intermediate regime, whose size is proportional to 
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FIG. 5: (color online) Result of the box counting algorithm 
for N = 262144, t « 45fo and varying volume fraction 
= 15.6%, 7.81%, 3.90%, 1.95%, 0.98%, 0.49% (from left to 
right); the lines along the data points are fits; the outer lines 
have slope —2 (top) and —3 (bottom); inset: scaling plot of 
Nhox/ffy^ versus Lb OX 

Results of the box counting algorithm for different den- 
sities (all in the percolating regime) are presented in 
Fig. [HI One can clearly distinguish two scaling regimes: 
at intermediate length scales a nontrivial slope of —2 
is observed, before the crossover to the expected value 
of —3 on the largest scales takes place. The size of 
the nontrivial regime shrinks with increasing density ac- 
cording to |log(/)|. Thus the crossover length scales as 
Lbox/c^ oc and the number of boxes at this length 
scales is A^box = (L/Lbox)'^ oc 0^. Hence if we plot 
A^box/0^ versus Lbox</'/rf, the data should collapse onto a 
single universal curve. This is indeed the case, as can be 
seen in the inset of Fig. [5l 

Long time asymptotics — In the late stage of cooling the 
probability Pbi, to break a bond is given by the probability 
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FIG. 6: (color online) Asymptotic time dependence for several 
volume fractions as in Fig.([l|; data (dots) in comparison to 
the analytical results (lines) 

to find a kinetic energy larger than AE: 

Pf,i,^ J d^ve{mv^/2~ AE)w{v). (5) 

We approximate the velocity distribution, w{v), by a 
Maxwellian of width T{t) and evaluate the above inte- 
gral in the hmit T{t)/AE 0. The probabihty to break 
a bond becomes exponentially small in that limit 

n.= (^) e--/- (6) 



dramatically slowing down the decay of the kinetic energy 
according to Eq.dl]): ^ ~ exp{—AE/T), predicting a 
logarithmic decay of T for asymptotically large times. 
The full solution of the above equation is in good agree- 
ment with the simulation data as shown in Fig. [6] The 
proportionality constant in the rate equation is fitted, be- 
cause the collision frequency is not known in this highly 
clustered state. 

Conclusions — The characteristic energy loss AE of wet 
granulates gives rise to a well defined dynamic tran- 
sition at time when the kinetic energy T is equal 
to AE. For t < to very few collisions give rise to 
bound pairs of particles and cooling proceeds according 
to T{t) = To (1 - t/to)"^. For t > to the particles aggre- 
gate into clusters with a flocculation regime in the early 
stages of aggregation and subsequent percolation. Clus- 
ter growth is a self similar process with a cluster size dis- 
tribution that obeys scaling. The asymptotic time decay 
is logarithmically slow, because bond-breaking becomes 
a very unlikely process. 
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